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ABSTRACT 

The observed properties of supermassive black holes suggest a fundamental link between their assembly 
and the formation of their host spheroids. We model the growth and activity of black holes in galaxies 
using ACDM cosmological hydrodynamic simulations by following the evolution of the baryonic mass 
component in galaxy potential wells. We find that the observed steep relation between black hole 
mass and spheroid velocity dispersion, Mbh cx a 4 , is reproduced if the gas mass in bulges is linearly 
proportional to the black hole mass. To a good approximation, this is equivalent to assuming the 
conversion of a fixed fraction of gas mass into black hole mass. In this model, star formation and 
supernova feedback in the gas are sufficient for regulating and limiting the growth of the central black 
hole and of its gas supply. Black hole growth saturates because of the competition with star-formation 
and in particular feedback, both of which determine the gas fraction available for accretion. Unless 
other processes also operate, we predict that the Mbh — cr relation is not set in primordial structures 
but is fully established at low redshifts, z < 2, and is shallower at earlier times. Once this relation is 
established we find that that central black hole masses are related to their dark matter halos simply via 

Mbh 9S We assume that galaxies undergo a quasar phase with a typical lifetime, £q ~ 2 x 10 7 

yr, the only free parameter of the model, and show that star-formation regulated depletion of gas in 
spheroids is sufficient to explain, for the most part, the decrease of the quasar population at redshift 
z < 3 in the optical blue band. However, with the simplest assumption of a redshift independent quasar 
lifetime, the model slightly overpredicts optical quasar numbers at high redshifts although it yields the 
observed evolution of number density of X-ray faint quasars over the whole redshift range 1 < z < 6. 
Finally, we find that the majority of black hole mass is assembled in galaxies by z ~ 3 and that the 
black hole accretion rate density peaks in rough correspondence to the star formation rate density at 
z ~ 4 - 5. 

Subject headings: accretion — black hole physics — galaxy: evolution — methods: numerical 



1. INTRODUCTION 

Recent dynamical studies indicate that supermassive 
black holes reside at the centers of most nearby galax- 
ies (e.g., Kormendy & Richstone 1995; Richstone et al. 
1998; Kormendy & Gebhardt 2001). The evidence indi- 
cates that the mass of the central black hole is correlated 
with the bulge luminosity (e.g., Magorrian et al. 1998) 
and even more tightly with the velocity dispersion of its 
host bulge, where it is found that A^bh oc a 4 (Tremaine 
et al. 2002; Ferrarese & Merritt 2000; Merritt & Ferrarese 
2001; Gebhardt et al. 2000). The tight relation between 
the mass of the black holes and the gravitational potential 
well that hosts them strengthens the theoretical arguments 
that there is a fundamental link between the assembly of 
black holes and the formation of spheroids in galaxy halos. 
In addition, the locally inferred black hole mass density ap- 
pears to be broadly consistent with the density accreted 
during the quasar phase (e.g; Soltan 1982; Fabian & Iwa- 
sawa 1999; Yu & Tremaine 2002) further supporting the 
idea that the nuclear activity, the growth of the black holes 
and spheroid formation are all closely linked. 

It is not yet determined what fundamental physical pro- 
cess is responsible for establishing the observed Mbh — cr 
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relation or what drives the strong evolution of quasars at 
low redshifts. Suggestions for the Mbh — a dependence 
have appealed to strong feedback due to quasar outflows on 
the protogalactic gas reservoir (e.g., Silk & Rees 1998); to 
capture of stars following protogalactic collapse (Adams, 
Graff & Richstone 2001) and to star formation regulated 
black hole growth (Burkert & Silk 2001). Models within 
the context of cold dark matter (CDM) cosmologies, in 
which the black hole growth is linked to the growth of 
galactic halos and activity is triggered during major merg- 
ers often make use of the observed Mbh — c relation (e.g., 
Wyithe & Loeb 2002; Volonteri, Haardt & Madau 2002 
and references therein) , or derive it dependent on some ad- 
ditional model assumptions (e.g., Kauffmann & Haehnelt 
2000; Haehnelt & Kauffmann 2001). 

In the commonly adopted merging scenario for the for- 
mation of spheroids, it is well established that black hole 
masses resulting from the coalescence of progenitor black 
holes fall short of the measured values (e.g., Ciotti & van 
Albada 2001). A large fraction of the black hole mass 
must therefore have been accreted, presumably as gas. At 
the same time, the gas present in the central region of 
spheroids must provide the gas supply that will eventually 
be accreted onto the black hole. In the simplest scenario 
the growth of the central black hole and the amount of fuel 
available (hence the black hole activity) should therefore 
be tied to the evolution of the gas component in galaxies. 

In this paper we investigate the growth and activity 
of black holes in bulges and attempt to look for correla- 
tions with the large scale properties of galaxies that match 
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those of the observed black holes at their centers. We use 
ACDM cosmological hydrodynamic simulations to follow 
the evolution of the baryonic mass in the potential wells of 
galaxies. The simulations (Springel & Hernquist 2003a, b) 
include a new prescription for star formation and feed- 
back processes within the interstellar medium which has 
been shown to produce a numerically converged predic- 
tion for the cosmic star formation history. Our aim is to 
trace from early times the gas fraction available for accre- 
tion and hence for the growth and fueling of black holes 
in competition with star formation and spheroid growth. 
This should allow us to determine directly to what extent 
both the Mbh — o relation for black holes, and the de- 
cline of quasars at low redshifts, are established through 
self-regulated star formation processes in the popular cold 
dark matter cosmogonies. 

In §2 we briefly describe the simulations and our anal- 
ysis. In §3 and §4 we present our results for the Mbh — c 
relation and our predictions for the luminosity function 
and space density evolution of quasars. We discuss fur- 
ther implications of our results in §5. In this paper we 
will limit ourselves to redshifts z < 7 where the quasar 
luminosity function has been measured. In future work 
we plan to extend our analysis to make predictions for the 
QSO population up to very high redshifts. 

2. SIMULATIONS AND ANALYSIS 

Throughout, we shall use a set of cosmological simu- 
lations for a ACDM model, with = 0.7, £! m = 0.3, 
baryon density fib = 0.04, a Hubble constant H = 
lOO/ikms -1 Mpc -1 (with h = 0.7) and a scale-invariant 
primordial power spectrum with index n = 1, normalized 
to the abundance of rich galaxy clusters at the present 
day (erg = 0.9). We here briefly summarize the main fea- 
tures of our simulation methodology and refer to Springel 
& Hernquist (2003a, b) for a more detailed description. 

Besides self-gravity of baryons and collisionlcss dark 
matter, the simulations follow hydrodynamical shocks, 
and include radiative heating and cooling processes of a 
primordial mix of helium and hydrogen, subject to a spa- 
tially uniform, time-dependent UV background (see, e.g., 
Katz ct al. 1996, Dave et al. 1999). The dark matter and 
gas are both represented computationally by particles. In 
the case of the gas, we use a smoothed particle hydrody- 
namics (SPH) treatment (e.g., Springel et al. 2001) in its 
entropy formulation (Springel & Hernquist 2002). 

In our numerical model, we assume that collapsed gas 
at very high overdensity becomes available for star forma- 
tion, which in turn creates a complex multi-phase inter- 
stellar medium (ISM). We use a sub-resolution model to 
describe the structure and dynamics of the ISM on unre- 
solved scales. In this model, star formation is assumed to 
occur in cold clouds that form by thermal instability out of 
a hot ambient medium, which is heated by supernova ex- 
plosions. These supernovae also evaporate clouds, thereby 
establishing a tight self-regulation cycle for star formation 
in the ISM. 

Motivated by observational evidence for the ubiqui- 
tous presence of strong galactic outflows in actively star- 
forming galaxies (e.g.; Martin et al. 1999; Heckman et al. 
2000) and by the required enrichment of the the interstel- 
lar medium (together with the relatively low efficiency of 
the global star formation) we have also included a phe- 



nomcnological treatment of galactic winds. In the fidu- 
cial parameterization of this process, each star forming re- 
gion is assumed to drive a wind with a mass outflow rate 
twice the star formation rate. We then assume that the 
kinetic energy of the galactic winds is comparable to the 
total available energy released by the supernovae associ- 
ated with star formation. This parametrization leads to an 
initial wind speed in the simulations equal to 484 km s -1 . 
These wind parameters are chosen to be representative 
of typical properties of the outflows associated with star 
forming disks (e.g. Martin, 1999). We note that whether 
or not the wind will escape from a galaxy (as it interacts 
and entrains infalling gas, shocks the ISM etc.) depends 
primarily on the depth of the galaxy potential well. For 
our choice of wind parameters we expect that only those 
halos with virial temperatures below 10 6 K have central 
escape velocities lower than the wind speed and will then 
lose some baryons in the outflows (c.f. §3.2). 

Note that our simulations did not include any prescrip- 
tion for black hole growth when they were run. The gas 
in our simulated galaxies is therefore only affected by star 
formation and winds. 

In order to resolve the full history of cosmic star forma- 
tion, Springel & Hernquist (2003a, b) simulated a series of 
cosmological volumes with sizes ranging from lh^ 1 Mpc 
to 100/i -1 Mpc. For each box size, a series of simulations 
was computed where the mass resolution was increased 
systematically in steps of 1.5 3 , and the spatial resolution 
in steps of 1.5. 

In this work, we show results mostly from their 'D- 
Series' which employed a periodic box of comoving side 
length 33.75 /i" 1 Mpc. Within this series, we use the 'D3', 
'D4' and 'D5' simulation runs with resolutions of 2 x 144 3 , 
2 x 216 3 and 2 x 324 3 particles, respectively, correspond- 
ing to mass resolutions in the gas of 1.43 x 10 8 /i~ 1 M Q , 
4.24 x lO 7 /^ 1 M and 1.26 x lO 7 /^ 1 M . The spatial res- 
olution of these simulations can be characterized by their 
gravitational softening lengths, which are equal to 9.38, 
6.25, and 4.17 h~ 1 kpc (in comoving units), respectively. 
The simulations of the D-Series have only been evolved to 
a minimum redshift of z = 1, because at this epoch, their 
fundamental mode starts to become non-linear. Among 
the simulation program studied by Springel & Hernquist 
(2003b), the D-Series constitutes the best compromise be- 
tween redshift coverage and resolution appropriate for the 
work here. The analysis of three simulations of different 
resolution but equal volume allows us to cleanly assess 
numerical convergence for all the physical quantities we 
measure, which is a significant advantage of our method- 
ology. For example, Springel & Hernquist (2003b) show 
that the star formation has fully converged below z ~ 6 
in the D5 run, while this is only the case for z ~ 5 in 
D4, and for z ~ 4 in D3. In general, we therefore expect 
that progressively higher resolution will also be required 
to achieve convergence for measurements of other physical 
quantities related to dissipative processes in the gas. We 
hence analyze simulation outputs for a set of 10 redshifts, 
given by z = 7, 6, 5, 4, 3.5, 3, 2.5, 2, 1.5, and 1. 

Note that in this work we attempt to link the large-scale 
properties of galaxies to those of observed black holes; for 
this goal it is not important to numerically resolve the 
properties of the accretion flows close to black holes. 
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Fig. 1. — The evolution from z = 7 to z = 1 of total gas mass, M gas , versus circular velocity V c , for galaxies in the D5 simulation. The 
dashed line shows the mean relation. The solid line shows a power-law (in log) fit to the points. The power-law index of the fitted line varies 
from 2.95 at z = 7 to 3.98 at z = 1 (sec text). 



2.1. Galaxy definition 

For the galaxy selection we used the Friends-of-Friends 
(FoF) algorithm to find groups of particles in the various 
simulation outputs. A linking length of 0.1 times the mean 
spacing of dark matter particles was used (also equal to 
0.1 x 2 1 / 3 = 0.125 times the mean spacing of the initial to- 
tal number of particles). The algorithm was applied to all 
star, gas and dark matter particles. In choosing the linking 
length, we have ignored the increase in the total number of 
particles due to creation of star particles as the simulation 
evolves. Making allowances for this would change the link- 
ing length by about 1.5% in the most extreme case. We 
discarded groups that contained fewer than 32 particles 
and also those with either no gas or no star particles. The 
FoF algorithm was chosen in the interests of simplicity, 
but we have checked that using another groupfinder such 
as SKID (e.g., Governato et al. 1997) makes no significant 
difference to our results. 

For each redshift output, we obtain a list of halos with 
their stellar (M star ), baryonic (M h — M gas + M star ) and 
dark matter (Mdm) mass components, and their positions. 
The total mass (M = Mb + Mdm) of each galaxy is used to 
give a measure of its radius and circular velocity. Following 
Mo & White (2002) and Springel & Hernquist (2003b), we 
assign a physical radius to a group of mass M according 
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and a corresponding circular velocity V c — (GM/r 5 oo) 1 ^ 2 , 
which may be expressed as 
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In these equations, we assumed that the selected groups 
on average enclose an overdensity of 500 times the back- 
ground value, which is slightly higher than the canonical 
value of 200 usually adopted in dark matter only simula- 
tions when a linking length of 0.2 is employed. In this way, 
we compensate for the bias in the measured group masses 
that arises from the slightly different distributions of gas 
and dark matter particles on the scale of the 'virial radius' 
^200 • Within this radius, the gas particles are actually 
slightly underabundant compared to dark matter parti- 
cles, such that we effectively select groups at a somewhat 
higher overdensity contour compared to a corresponding 
dark matter only simulation. In fact, we have directly 
compared our measured group masses to an alternative 
selection that defines groups based on FoF applied to the 
dark matter only, using the canonical linking length of 
0.2, and increasing the final group mass by the enclosed 
gas and star particles. In this case, the density contrast 
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Fig. 2. — The 3D velocity dispersion of star particles, a, versus 
circular velocity, V c , of galaxies in the D5 simulation output at z = 
1. The solid line is a linear fit (in log space) to the points. 

with respect to the background can be assumed to be 200, 
and the group mass M200 may be defined to represent the 
'virial mass'. We have correlated the two mass definitions 
using an object-to-object comparison and found very small 
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Fig. 3. — The relation we obtain between Mbh (assumed to be 
proportional to the total gas mass) and velocity dispersion a of stars 
for groups in the 'D5' simulation. The dashed line shows the mean 
relation over all redshifts, i.e. for z = 7, 6, 5, 4, 3.5, 3, 2.5, 2, 1.5, 
and 1. For comparison, the solid line shows a power-law with a 
logarithmic slope of 4. 



scatter between the two different selection methods, with 
the two masses being related by 

M 200 ^ 1.6M, (3) 

which can be used to approximately convert the masses we 
quote in this paper to the virial masses of halos. We will 
also use the M 2 oo catalogues to compare the differences 
between the baryon fraction in objects selected with the 
two different methods (see §3.2). 

For each galaxy we also directly calculate 3D stellar ve- 
locity dispersions, a, from the star particles. Finally, we 
measure the amount of cold gas M co id in galaxies, which 
we here define to be the mass of actively star-forming cold 
clouds in the ISM. 

3. COSMOLOGICAL BLACK HOLE GROWTH 

In Figure 1, we plot the gas mass M gas for each in- 
dividual object selected in the six simulation outputs at 
redshifts z = 7, 5, 4 , 3, 2, and 1, versus its circular ve- 
locity V c . The dashed line shows the average relationship 
traced by the galaxies at the different redshifts and the 
solid line represents the best power law fit to the objects 
with V c > 80kms _1 . The fit has functional form 

where the constants a and b for the six redshifts are (z = 
7; 2.95, -9.02); (z = 5; 3.00, -9.02); (z = 4; 3.08, -9.12); 
(z = 3; 3.10, -8.98); (z = 2; 3.65, -10.2); (z = 1; 3.98, -10.4) 
Here we do not attempt to assess the significance of the 
correlations as the sources of scatter in the measurement 
of M gas cannot be easily quantified (see §3.1); the fitting 
is merely intended to illustrate the overall evolution of 
the slope, a, that best describes the relation. The most 
remarkable feature that we note is the steepening of the 
-Mgas — V c relation from values a ~ 3 at high redshifts, to 
a ~ 4 at z < 2. 

We then compare the slopes (in log) of the M gas — V c 
relation to the most recently revised Mbh — & relation from 
Tremaine et al. (2002) 

/ a \ 4.02±0.32 

M BH = (1.5 ± 0.2) x 10 8 M Q ( r ) . (5) 

V ; \200kms- 1 J V ' 

(Note that the first published estimates of the slope (in 
log) of the above relation varied from values ~ 4.8; Fer- 
rarese & Merritt (2000); Merritt & Ferrarese (2001); to 
- 3.75 Gebhardt et al. (2000), highlighting the difficulty 
of these measurements). In our lowest redshift output, at 
z = 1, the observed slope of the Mbh — o relation (Eq. 5) 
is reproduced well by the M gas — V c relation (Equation 4), 
suggesting that the slope of Equation 5 is easily explained 
if Mbh is directly proportional to M gas (similar slopes are 
obtained if we restrict ourselves to V c < 400 km s -1 , con- 
sistent with the observed velocity range used to derive the 
Tremaine et al. 2002 relation). At increasingly higher 
redshifts, the M gas — V c relation flattens and approaches 
M gas oc V c 3 at z ~ 7, as shown by the fitted slopes in 
Equation (4) and as expected from Equation (2). 

To make more direct contact with the observed Mbh — 0" 
relation and to consider its possible link with the M gas — V c 
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relation in the simulations, we need to consider the rela- 
tionship between V c and the stellar velocity dispersion a, 
as measured directly from the simulations. The results 
for the a — V c relationship measured in the simulation are 
shown in Figure 2; the plot shows the z = 1 simulation out- 
put, although we note that we find no significant temporal 
evolution in this relation. There is a direct proportionality 
between V c and a,^ which we fit with the linear relation 
(in log): 

log cr = 1.03 log 14-0.32, (6) 

shown by the solid line in Figure 2. Because the mea- 
surement of a from the star particles is intrinsically much 
noisier than that of V c (which is effectively a measure of 
mass), we make use of the relation from (Eq.6) to con- 
vert from V c to a in what follows. Using a directly would 
introduce a significant source of artificial scatter in our 
simulated Mbh — o relation. We note that in future work 
it might be possible to investigate the scatter about the 
Mbh — a relation provided it is possible to eliminate nu- 
merical scatter in the simulated stellar velocity dispersion 
a. 

The two relations, M gas — V c and V c — a from the simula- 
tions, if compared to the observed Mbh — cr relation (Eq. 5) 
suggest that the central black hole mass in galaxy may be 
simply proportional to the gas mass in galaxies. We hence 
make this simple ansatz that the black hole mass is pro- 
portional to the gas mass, allowing us to plot a predicted 
relation between Mbh and a obtained from the simula- 
tions. In Figure 3, we overplot the corresponding data 
from all the redshift outputs used in Figure 1, resulting in 
a composite average relation which is shown as a dashed 
line. Only simulation outputs for the D5 run are included, 
and we defer a discussion of differences between the D4 
and D5 runs to the next section. In order to match the 
normalization of the observed Mbh — cr relation (Eq. 5) we 
associate black hole and gas masses here according to 



M BH = / M g 



0.004/i" 1 M e 
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where the solid line in Figure 3 shows Equation 5. Assum- 
ing that there is no significant change from z = 1 to z = 0, 
(which we have checked is a good approximation by mak- 
ing use of the larger 100 h _1 Mpc simulation box in the 
G-Series; see Springel & Hcrnquist 2003a) this therefore 
implies that in the context of our model, the mass of the 
central black hole in galaxies is consistent with being < 1 
per cent of the total gas mass in galaxies. 

In this interpretation, the Mbh — cr relation is not set 
in primordial structures and maintained through cosmic 
time; it is established at later times as the growth of 
black holes is limited by the available gas, in competition 
with star formation and associated feedback processes like 
galactic winds. Because both the black hole growth and 
the change of the amount of available gas occur gradually 
from high to low rcdshifts, the resulting correlation does 
not show large scatter. Note, however, that in this simple 
model, the high-mass, high-redshift objects are predicted 
to lie on the right of the main relation. Large samples 

tThe direct proportionality between V c and cr also suggests that 
our group definition methods described in §2.1 should be a fairly 
reliable in their identification of galaxies. Note that the scattering 
of points above the mean relation in Figure 2 is most likely due to 
merging systems. 



of intermediate mass objects at z > 3 will be required to 
determine observationally whether the Mbh — cr relation 
is indeed shallower at earlier times, as predicted here. At 
present the observed Mbh — cr relation (as shown in Eq. 5) 
is derived from local objects. For objects at z = 1 (sixth 
panel of Fig. 1) the predicted relation is extremely tight. 

Note also that given the cr — V c relation (Eqs. 6 and 7) 
from the simulations we find that the black hole mass is 
related to the dark matter mass of the host galaxy by the 
relation: 



M BH 
10 8 M« 



0.7 



M DM 



4/3 



(8) 



This is broadly consistent with the Mbh — V c relation de- 
rived by Ferrarese (2002) and Baes et al. (2003), reflecting 
the fact that the relation in Equation 6 matches the ob- 
served correlation between V c and a (Ferrarese 2002; Baes 
et al. 2003). 

3.1. Numerical convergence and scatter 
3.1.1. Deviations from a power law 

Figures 1 and 3 show that, for small values of M gas 
and V c , the Mbh — cr relation deviates significantly from 
a power-law and its scatter is increased. Figure 4 shows 
the same plot as in Figure 1 but with the D3, D4 and 
D5 simulations over-plotted at z = 1 in the top panel, 
and the D5 and Q5 simulations in the bottom panel. The 
'Q5' simulation is taken from the Q-Series of Springel & 
Hcrnquist (2003b). While Q5 has a smaller cosmological 
volume of size 10/i -1 Mpc on a side, its particle number is 
equal to that of D5, so that the mass resolution in the gas 
becomes 3.26 x 10 5 /i _1 M Q (a factor of 40 improvement 
over D5), and the spatial resolution becomes 0.63 hr 1 kpc 
in comoving units. 

In the top panel of Figure 4, the average relations for 
the D3, D4 and D5 runs are shown with the dashed, solid 
and dot-dashed lines, respectively. The figure shows that 
both the scatter and the deviation from the relation de- 
crease with increasing resolution of the simulation. This 
is a result of the fact that close to the resolution limit of 
each of the simulations, star formation cannot be resolved 
properly any more, making the measurement of M gas un- 
reliable. For higher resolution simulations, this limit for 
convergence in M gas moves to lower mass systems. 

However, when analyzed in isolation, even for the D5 
run it is not clear whether the flattening of the relation 
at small V c is due to some underlying physical process as- 
sociated with star formation, or whether it can be fully 
explained by insufficient resolution. We investigate this 
further in the bottom panel of Figure 4, where we show a 
comparison of the measured gas content between the D5 
and the Q5 simulations at redshift z = 3 (this is close to 
the minimum redshift reached by the Q-Series). Although 
the number of objects in the high mass range is much 
reduced in Q5 due its smaller volume, the increased reso- 
lution allows us to investigate the M gas — V c relation for 
much smaller objects. We find almost negligible deviations 
from the extrapolated power law behavior derived from the 
D-Series. It is therefore plausible that the Mbh — M gas re- 
lationship extends to very low V c given sufficient numerical 
resolution. 

Determining the extent to which our theoretical relation 
applies to the low mass range is also relevant in view of re- 
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Fig. 4. — Convergence study for the measurement of M gas . We 
plot the same quantities as in Figure 1. Top panel: groups extracted 
from the D3, D4 and D5 simulations at rcdshift z = 1 arc overplot- 
ted with plus signs, triangles and cross symbols, and in shades from 
dark to light gray, respectively. The dashed line is the average re- 
lation in the D3 run, the solid line in D4 and the dash-dotted in 
D5. The large scatter at small V c is reduced with increased reso- 
lution. Bottom panel: We compare here objects in the Q5 (dashed 
line) and D5 simulations at rcdshift 2 = 3. Thanks to the sig- 
nificantly increased resolution of the Q5 simulation, the scatter is 
further reduced, showing that it largely arises from resolution ef- 
fects. However, the small box of Q5 contains only a few high-mass 
objects. 



cent searches for intermediate mass black holes in globular 
clusters. The promising detection of a central black hole 
of M BH ~ 2 x 10 3 M q in stellar cluster Gl in M31 (Geb- 
hardt, Rich & Ho, 2002) and the ongoing debate on the 
presence of a similar central object in the globular cluster 



M15 (Gerssen et al. 2002) seem to be consistent with an 
extrapolation of the observed Mbh — f relation for galax- 
ies, possibly implying a similar formation process for the 
central black hole. 

Figures 1 and 3 also show that for large values of V c the 
relationship between M gas and V c flattens and almost re- 
covers the power law slope of ~ 3 expected from Equation 
(2). for halos with baryon content consistent with the uni- 
versal fraction. As we will discuss further in §3.2, stellar 
feedback with its associated outflows plays a fundamental 
role in our model for generating the steep decline of gas 
content towards low-mass halos, but these processes even- 
tually become inefficient for very massive systems, which 
then recover the universal baryon content. 

3.2. The evolution of M gas , M star and M co \^ 

In order to establish the constant of proportionality 
/ in Equation (7) we need to investigate whether the 
measurement of M gas has converged for sufficiently large 
V c > 100 km s -1 . Figure 5 shows the evolution of the to- 
tal gas, stellar and cold gas (gas subject to star formation) 
masses, M gas , M sta r and M co id as a function of z at con- 
stant circular velocity, V c = 200kms _1 . We also include 
the dark matter component, Mdm, for comparison. The 
plot hence shows the evolution of the baryonic (and dark 
matter) content in typical L*-galaxies. The black symbols 
and lines are for the D3 simulation, the light gray for the 
D4 run and the dark gray for the D5 run. For each simula- 
tion, the star symbols joined by dot-dashed lines represent 
M star , the triangles joined by the dashed lines M col d, and 
the solid dots joined by the solid lines M gas . Mdm is shown 
for the D5 run only, using diamonds connected by a dotted 
line and rescaled by a factor £lb/Qr>M = 0.153. 



The figure shows that the measurements of M g 



M s1 



and M co id at V c = 200 km s -1 have converged in the D4 
run to good accuracy. In the D3 simulation, the star for- 
mation has not fully converged even at this circular veloc- 
ity; hence M star is lower in the D3 run than in D4 or D5, 
and M co id, in particular, as well as M gas are overestimated. 

It is further seen that the ratio M gas /M star is large at 
high redshifts and approaches values > 1 at redshifts z < 
3. The cold gas fraction M co id/M star is also fairly large at 
high rcdshift but decreases significantly below 1 at z ~ 1. 
Interestingly, in this plot we see the overall evolution of 
the baryon and dark matter mass in galaxies. This will 
be particularly useful when deriving the evolution of the 
quasar luminosity function in §4. Here we find that M gas 
for the circular velocity given, V c = 200kms _1 , is a slowly 
decreasing function of redshift, where we have M gas oc (l + 
z)~ 1 / 3 , as opposed to Mdm, which is a steeper function 
of redshift, following the relation Mdm oc (1 + z)~ 3 / 2 as 
expected from Equation (2). Even though gas is both used 
up by star formation and expelled by galactic outflows, the 
total mass of gas does not decrease with decreasing redshift 
but remains nearly constant. 

Note that in the Kauffman & Haehnelt (1999) model, in 
order to explain the small scatter in the Mbh - o relation, 
Haehnelt & Kauffmann (2000) require the cold gas mass 
of their bulge progenitors to be essentially independent 
of redshift (where in their model a fixed fraction of the 
cold gas is assumed to to be accreted by the black holes). 
Here we find that M co id decreases with redshift and instead 
the total gas mass M gas is found to be a slowly varying 
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Fig. 5. — Gas, star, cold gas and dark matter mass components 
for halos of constant circular velocity Vc = 200 km s -1 as a function 
of redshift. The results for the D3 run are shown in black, for D4 in 
light gray, and for D5 in dark gray. The dark matter mass has been 
rescaled by a factor f2(,/f2DM = 0.153 and is shown for comparison 
only (for D5). Approximate convergence of the various baryonic 
mass components is already obtained for the resolution of the D4 
simulation. 



function of redshift. This suggests that the slope of the 
-Mbh - cr relationship may not be set by the amount of 
gas that cools and becomes available for star formation, 
as assumed by Haehnelt & Kauffmann (2000), but rather 
by the total amount of gas in the galaxies, although it is 
certainly true that galaxies have higher cold gas fractions 
at earlier times. In the next section we investigate how 
star formation and feedback regulate the total amount of 
gas available for accretion onto central black holes. 

3.3. The role of stellar feedback 

In Figure 6, we plot the baryonic mass fraction (M gas + 
M stax )/M tot in galaxies as a function of V c for the objects 
in the D5 simulation at the three redshifts z = 5, 3, and 1, 
overplotted in black, dark gray and light gray, respectively. 
A running average of the baryon fraction at each redshift 
is shown as a dash-dotted line for z = 5, dashed line for 
z = 3 and solid line for z = 1. The solid horizontal line in- 
dicates the universal baryon fraction f2f,/f2 m expected for 
our cosmological model. The two panels of the figure show 
a comparison of the baryon fractions measured for objects 
selected according to the two group definitions outlined in 
§2.1. The top panel uses the same definition as in our pre- 
vious figures, where we applied the FoF algorithm to dark 
matter, gas and star particles on an equal footing, while in 
the bottom panel, we selected objects based on dark mat- 
ter particles only, and then included all enclosed gas and 
star particles. Based on this comparison, we note that the 
larger scatter in the top panel at low V c is a consequence of 
the different group selection methods. In particular, close 
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Fig. 6. — The baryon mass fraction in galaxies as a function of 
circular velocity V c . Three redshifts, z = 5, 3, and 1, are shown 
in shades of dark to light gray, and we plot running averages for 
them in dot-dashed, dashed and solid lines, respectively. The hor- 
izontal line shows the baryon fraction f!b/fi m . In the top panel, 
objects were selected by applying a FoF linking length of 0.1 the 
initial dark matter particle spacing to all particles, with the total 
mass in galaxies defined as Mtot = M gas + M s tars + JWdm- ln the 
bottom panel, a linking length of 0.2 has been used instead, but was 
restricted to the set of dark matter particles. Here, the enclosed gas 
and star particles in each identified group have been added in a 
separate step (see § 2.1). 



to the resolution limit our default method is somewhat 
more prone to occasionally producing measured baryon 
fractions well above Q(,/fi m , as it can happen if the FoF 
algorithm just picks out the concentrated cold gas and 
stellar parts of low-mass galaxies, but loses the sparsely 
sampled dark halo around them. Although this occurs, 
our method is able to recover the universal baryon frac- 
tion at z = 1 in the simulation more closely (top panel, 
Fig. 6) than in the bottom panel. 
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However, more importantly, both plots show that the to- 
tal baryon mass fraction decreases from intermediate val- 
ues of V c towards low values, while systems of very large 
mass remain close to the universal baryon fraction. Fur- 
thermore, the strength of the baryonic depletion at low V c 
increases towards low redshifts. 

It is clear that these trends in Figure 6 are a reflection 
of the same steeping of the M gas — V c relation towards low 
redshift and low V c that we already observed in Figure 1. 
However, in this representation the fundamental physical 
mechanism that drives this evolution is more apparent. In 
essence, the M gas — V c relation is steepened by the loss of 
baryons due to strong feedback by galactic winds. These 
galactic outflows are increasingly more efficient in low mass 
systems (where resolved; see also Springel & Hernquist 
2003b), while for very massive groups, the winds eventu- 
ally become trapped in deep gravitational potential wells. 
Note that the measured M gas -V c relation indeed flattens 
for high-mass systems (as shown by Figures 1 and 3), with 
a transition scale that is consistent with the maximum cir- 
cular velocity expected to allow an escape of winds from 
halos in the simulations. 

4. BLACK HOLE ACTIVITY 

Quasars have long been believed to be powered by the 
accretion of gas onto supermassive black holes and have 
been known to evolve very strongly. Their comoving space 
density increases by nearly two orders of magnitude out to 
a peak at z <~ 3 (e.g., Boyle et al. 2000). Their evolution 
at higher redshifts is more uncertain although recent ob- 
servations have begun to constrain the bright end of the 
luminosity function out to z ~ 5 and individual quasars 
have been detected out to z — 6.28 (Fan et al. 2001a,b). 
Here we use our results from §3 to predict the evolution of 
the number density of bright QSOs from redshift z = 7 to 
z = 1 and compare it with observations. 

4.1. Gas accretion history 

In order to predict black hole luminosities, we need to 
estimate the amount of gas available for accretion for a 
given characteristic accretion timescale. In §3 we showed 
that the mass of the black hole may simply be modeled 
as a fixed fraction of the total amount of gas in the halo, 
suggesting that black holes grow primarily by accretion 
of gas, in which case the amount of accreted mass also 
depends on star formation and feedback processes. 

The small-scale physical processes linked to the growth 
of black holes, however, do not give us a direct handle 
for constraining the black hole accretion timescale. Given 
the lack of any definite prescription or constraint for the 
quasar lifetime or its evolution, we make the simplest pos- 
sible assumption and choose the quasar lifetime which is 
independent of redshift and of order of the Salpeter time 
tq <~ 2 x 10 7 yr. Our choice of tQ is most appropriate 
for the bright end of the quasar population (for quasars 
that radiate around the Eddington limit) which is what we 
concentrate most of our modeling on, as the faint end is 
dominated by low mass galaxies which may not be fully re- 
solved in the simulations used here (see §3). We note that 
this choice of timescale is consistent with results of Steidel 
et al. (2002) who estimated the lifetime of bright quasar 
activity in a large sample of Lyman-break galaxies at z ~ 3 
to be ~ 10 7 yr and range of quasar lifetimes obtained 



from quasar clustering constraints ( which generally imply 
10 7 yr <t Q < 10 8 yr; e.g.; Martini & Weinberg 2001). It is 
interesting to note that a similar timescale, tQ ~ 2 x 10 7 
yr, was also adopted in the models of Sokasian, Abel and 
Hernquist (2002; 2003) developed to study reionization of 
Hell and the properties of the ionizing background at in- 
termediate redshifts. In this case, this timescale was shown 
to match the observations of the Hell Lyman alpha forest 
and of the H and Hell opacity at z ~ 2.5 — 5. 

As we shall see the choice of this parameter is crucial 
for determining the evolution of the quasar population. 
For our model, the star formation timescale is i star ~ a 
few Gyr; hence, we have tQ « t sta , T << Hq 1 , consistent 
with our hypothesis that the overall cosmological change 
in the environment governs the quasar activity and the 
duty cycle of quasars increases with redshift. 

We assume that all galaxies can undergo a quasar phase 
and so accrete a fraction of their total gas content. Be- 
cause, in general, the evolution of the gas mass in galaxies 
follows an evolution gS (1 + z) 1 / 3 (i.e.; 'on average' the gas 
mass does not go down, but merely does not grow; see 
Fig. 5), we can assume the black hole growth occurring in 
the last Hubble time always dominates. Hence, we con- 
sider the black hole mass at the time of accretion to be 
the contemporaneous black hole mass. This simple sce- 
nario is then consistent with the approximate condition, 
Mbh oc M gas relating the total gas and black hole masses. 

The bolometric luminosity of a quasar at redshift z is 
then simply given by 

L = V M(z)c 2 ~ v ^p c 2 ~ V fM saa (z)c 2 /t Q , (9) 

where we have taken AM gas (z) ~ fM gas (z) with / = 
0.4%, as suggested by Equation (7); i.e. the average mass 
accreted is a constant fraction of the total gas mass M gas . 
We adopt the standard value for accretion radiative effi- 
ciency of r\ = 10%, and assume that every object with a 
nonzero star formation rate (this discards small-mass ob- 
jects where star formation is likely not resolved) shines 
as a quasar for a time tQ. The expected number den- 
sity of active quasars at each redshift is a fraction /(zj) = 
tQ x (z i+1 - Zi)/[t H (zi) - t H (z i+1 )} where ^ and t H (zi) 
are the redshift and Hubble time of simulation output i, 
respectively. 

4.2. The evolution of the quasar luminosity function 

4.2.1. Optical blue band 

We now compare our simple predictions for the luminos- 
ity function to observational data. We take a fixed fraction 
eb of the bolometric luminosity to be radiated in the blue 
band. In the B-band, the bolometric correction cb, de- 
fined as L = cbvLb, is about 11.8 (Elvis et al. 1994), 
where vL b is the energy radiated at the central frequency 
of the -B-band per unit time per logarithmic interval of 
frequency At the peak of its lightcurve the quasar has a 
.B-band magnitude 

M B = 37.05 - 2.5 log(e Bj L), (10) 

where L is in units of ergs -1 . We place an upper limit 
on the maximum allowed luminosity of a quasar so that 
L does not exceed the Eddington luminosity L^dd- Note 
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Fig. 7. — The f?-band luminosity function of quasars at various redshifts, as indicated in each panel. The filled circles are the values from 
our models. Error bars indicate the Poissonian error in the counts. The solid lines and light gray shaded areas at z = 2, 1.5, and 1 show the 
best fit to the 2dF luminosity function (Boyle et al. 2000). The dark shaded areas and solid lines at z = 6, 4, and 3 give the best fit to the 
SDSS luminosity function from Fan et al. (2001a,b). 



that if accreted gas provides the power, then clearly there 
cannot be steady fueling at rates that generate L > L^dd- 
With the black hole masses given by Equation (7) a small 
percentage of objects radiate above -£<Edd- 

Figure 7 shows a comparison of our theoretical luminos- 
ity functions (LFs) at different redshifts with the most re- 
cent determinations of the quasar B-band LF from the 2dF 
QSO Rcdshift survey (0.3 < z < 2.3; Boyle et al. 2000) 
and the Sloan Digital Sky Survey SDSS (3.6 < z < 6; Fan 
et al. 2001a). For the 2dF, the luminosity function for 
—26 < Mb < —23 is described by a double power law 
with a bright end slope of /?i = —3.4 and a faint end slope 
of (3 = -1.5 (Boyle et al. 2000; Pei 1995). The QSO lu- 
minosity function measured by Fan et al. (measured for 
-25 < M1450 < -28; where M 1450 is the absolute AB 
magnitude; see Eq (5) in Fan et al. 2001a for the relation 
between Mb and Mi4 50 ) gives a flatter bright-end slope 
(/3i = —2.5). For the luminosity function at z = 3, we 
plot both these functions, one extrapolated from higher 
redshifts, the other from lower redshifts. 

Our simple model reproduces the luminosity function of 
optically selected quasars in the redshift range 4 < z < 1.5 
reasonably well. The overall turnover of the space den- 
sity of quasar at low redshifts (z < 3) is also reproduced 
with some overestimate of the number of bright quasars 
at z = 1. The slope of the LF at redshifts higher than 
z > 4 matches the steep value found by Fan et al. (2001a). 



However, the model at this high redshifts predicts a larger 
numbers of quasars (both bright and faint) with respect to 
observations (and their extrapolation to faint magnitudes 
M B > -25). At redshifts z = 3 and 2, the bright end 
of the luminosity function flattens in accordance with the 
observations by Boyle et al. (2000) but remains too flat at 
z = 1, leading to an overestimate of the number of bright 
quasars at these redshifts. 

In the top panel of Figure 8, we show the evolution 
of the space density of bright quasars with luminosity 
Mb < —26 as a function of redshift. In this plot, we add 
as open symbols data from the Warren, Hewett, & Os- 
mer survey (1994; hereafter WHO) in the redshift range 
2 < z < 4.5, and from Hartwick & Shade (1990) in the 
range 0.1 < z < 3.3, respectively (also summarized in Pei 
1995). The dotted line is a fit to the Schmidt, Schneider & 
Gunn survey (1995; SSG). The Fan et al. (2001a,b) data 
is shown with solid dots, while our theoretical prediction is 
shown using solid symbols joined by a solid line. Although 
this measurement is quite noisy as a result of a relatively 
small simulation box, the overall evolution of the comov- 
ing number density of quasars is reproduced reasonably 
well by the model. Note that because of our limited box 
size no objects can be found massive enough to contribute 
to the luminosity function at Mb < — 26 for z > 5. As 
we shall discuss in §4.2.2 , at high redshifts, black hole 
accretion follows the growth of structure; to a good ap- 
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proximation our prediction at higher redshifts should thus 
follow the extrapolation from lower redshifts. Because of 
the small number of objects and the very sharp decrease in 
the quasar counts for objects brighter than —26 our results 
are very sensitive to the model parameter tQ. In order to 
get a more robust comparison we compare our results to 
recent constraints on quasar number counts derived from 
recent Chandra X-ray observations which probe quasars of 
more intermediate luminosity. 

4.2.2. X-ray band 

The bottom panel of Figure 8 shows the redshift evo- 
lution of the comoving number density of quasar in the 
2 — 8 keV band compared to the recent Chandra observa- 
tions (Cowie et al. 2003; Barger et al. 2003). We plot 
our model with solid dots joined by a dashed line and a 
solid line for two observed luminosity ranges 10 43 < L x < 
lO^ergs" 1 and 10 43 < L x < lO^ergs -1 respectively. 
We take i x ,2-8keV ~ 0.03-L where L is the bolomctric lu- 
minosity given in Equation 9. This relation seems to be 
roughly valid for all type of AGN (Elvis et al. 1994). For 
the higher luminosity range we take tQ — 2 x 10 7 yr (as in 
§4.2.1), whereas for the lower luminosity range the quasar 
timcscale is a factor of 2 smaller, tQ = 1 x 10 7 yr, as this 
seems to describe better the overall evolution of these ob- 
jects. Overall, we find the model to agree pretty well with 
the current constraints from the X-ray band if we allow 
a small, mass dependent variation of the typical accretion 
timescale. This again emphasizes the strong dependence 
of our predictions on this parameter. 

4.2.3. Dust, obscured AGN, low-radiative efficiency 

Observations indicate that the comoving number den- 
sity of quasars peaks at z ~ 3 and drops abruptly by 
a factor of 20 from redshift z ~ 2 to z ~ 5. Our pre- 
dicted number density lies closer to an extrapolation of the 
low redshift luminosity function (e.g., Boyle et al. 2000), 
implying a perhaps less rapid evolution at high redshifts. 
From the observational point of view, there is still some 
uncertainty concerning the recent high-redshift results as 
a result of a number of potentially very important system- 
atic biases (see e.g., Pei 1995; Fan et al. 2001a). Surveys at 
high redshifts still suffer from incompleteness, uncertain- 
ties in the K-corrections (when compared to low redshifts) 
and in particular from the possibility that a large number 
of quasars are not detectable in optical surveys as a result 
of dust extinction. 

The tendency to overestimate the number of bright 
quasars at z < 1 in our model may also be less severe 
given the recent Chandra results, which imply the presence 
of a substantial population of optically obscured luminous 
AGNs at low redshifts (e.g., Rosati et al. 2002; Barger et 
al. 2001). This would help to reduce the discrepancy. We 
also note that there is significant evidence from Chandra 
observations of nearby galactic nuclei indicating that the 
gas supply onto their central black holes far exceeds their 
luminosity outputs, implying that the radiative efficiency 
r/ of nearby supermassive black holes may be a lot lower 
than the canonical value of r\ — 0.1 we have taken (e.g., the 
case of M87; Di Matteo et al. 2003; our Galactic Center, 
Baganoff et al. 2002; Narayan 2002). 

In summary, we have shown that the evolution of the 
gas supply via star formation and feedback can, for the 
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Fig. 8. — Top: The evolution of the space density of quasars with 
Mb < —26. Our model is shown by the solid circles joined by the 
thick solid lines. The model is for tQ = 2 X 10 7 yr. The dotted line 
is a fit to the SSG survey in the redshift range 2.75 < z < 4.75. The 
open squares show a compilation from Pei (1995) derived from the 
WHO and HS surveys. The filled circles are the SDSS points derived 
by Fan et al. (2001a, b). The shaded area at low redshift is the best 
fit model from the 2dF survey by Boyle et al. (2000). Bottom: The 
evolution of space density of quasar in the 2 — 8 keV X-ray energy 
band. The data are from recent Chandra observations by Cowie et 
al. (2003); Barger et al. (2003). The open squares represent the 
quasar number density of sources with 10 43 < L x < 10 44 crgs _1 , 
and the open triangle the number density of sources with 10 44 < 
L x < 10 45 ergs — 1 . The dotted bars show upper limits. The solid 
and dashed lines show our model in these luminosity ranges using 
with tQ = 2 X 10 7 yr and tQ = 10 7 yr respectively. 

most part and without the introduction of any additional 
parameters, reproduce the turnover in the comoving num- 
ber density of quasars at z < 3. We obtain an improved 
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agreement with observations with respect to previous mod- 
els in which the evolution of the luminosity function is 
predicted within the context of hierarchical build-up of 
the dark halos through merger trees (e.g., Kauffmann & 
Haehnelt 2000; Cattaneo 2000; Volonteri et al. 2002; Hatz- 
iminaoglou et al. 2002). However other factors such as 
obscuration or change in radiative efficiency or other feed- 
back processes (e.g.; due to radio jets) may still play some 
role in producing the very sharp decrease in number counts 
at z < 1. 

5. THE BLACK HOLE ACCRETION HISTORY 

The evolution of the mass density of black holes pre- 
dicted from our model is shown in the top panel of Figure 
9. Extrapolating to redshift z ~ 0, we obtain pbh ~ 2— 4x 
10 5 M Mpc" 3 , which is within 20 - 30% of the value ob- 
tained by Merritt & Ferrarese (2001), and is roughly con- 
sistent with the result of p B H = (2.5±0.4) x 10 5 M Q Mpc~ 3 
derived by Yu & Tremaine (2002) using the velocity dis- 
persions of galaxies in the SDSS (also similar to the value 
of Salucci et al. 1999). Values inferred from the X-ray 
background lie typically somewhat higher, in the range 
of pbh = 6 - 16.8 x 10 5 Mq Mpc~ 3 (Fabian & Iwasawa 
1999; Elvis, Risaliti & Zamorani 2002). Although the most 
recent X-ray constraint from Chandra observations give 
Pbh ~ 2 - 4M©Mpc~ 3 (Cowie et al. 2003). Note that 
Figure 9 clearly shows that the majority of the black hole 
mass is already built up by z ~ 4. 

In the bottom panel of Figure 9, we plot the evolution 
of the comoving black hole mass accretion rate density 
according to our model. As suggested by the cumulative 
black hole mass function in the top panel, the peak of 
the black hole accretion rate per unit volume occurs at 
rcdshifts around 4-5. 

It is also interesting to compare the black hole accretion 
rate density (BHAR) directly with the cosmic star forma- 
tion rate density (SFR), as measured by Springel & Hern- 
quist (2003b) for the simulation set we analyze here. The 
dashed line in the bottom panel of Figure 9 shows an em- 
pirical double-exponential fit used by Springel & Hcrnquist 
(2003b) to approximately describe the simulation results 
for star formation, here rescaled by a factor 2 x 10~ 3 . 

Comparing this fit with our results for black hole growth, 
it is evident that the shapes of the BHAR and SFR histo- 
ries are quite similar in general, particularly at high red- 
shifts, but that the BHAR history more steeply declines 
at low redshifts. 

We may hence also try to use the double exponential 
function employed by Springel & Hcrnquist (2003b) to fit 
the BHAR. Using the ansatz 



P{z) 



6cxp[q(^ - z m )} 
b — a + acxp[6(z — z 



(11) 



we obtain a good fit for our measurements of the BHAR 
with parameter values a = 5/4, b = 3/2, z m = 4.8 and 
e = £bh = 3 x 10 _4 M Q yr~ 1 Mpc -3 . This fit is shown as 
the solid line in the bottom panel of Figure 9. 

Compared to the SFR history, the peak of the BHAR ap- 
pears to occur at slightly lower redshift, but the high-z evo- 
lution is nevertheless broadly similar, suggesting that the 
evolution of both rates is largely driven by the rapid grav- 
itational growth of the halo mass function in this regime. 
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Fig. 9. — Top panel: The evolution of the black hole mass density 
in our model. The majority of black hole mass is assembled by 
z ~ 4. Bottom panel: The evolution of the black hole accretion rate 
(BHAR) density in our model for black hole activity described in §4. 
The solid circles mark the measurements from the simulations. The 
solid line is an analytical fit to our results using the fitting function 
of equation (8) and the dashed line using the function of equation 
(9) where the numerator is tx x 4 instead of \ 2 . The dashed-dotted 
and dotted lines show the fits to the cosmic star formation history 
in the simulation set analyzed by Springel & Hcrnquist (2003b) and 
Hcrnquist & Springel (2003), but scaled down by a factor of 500. 
In our model, the black hole accretion rate density has a peak at 
similar rcdshifts as the star formation rate density, but it decreases 
faster at low rcdshifts. 



Note however that even at its peak, the black hole accre- 
tion rate density is roughly a factor 500 lower than the 
star formation density, implying that the rate at which 
baryons are locked up in stars always far exceeds that in 
black holes. 
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It is also interesting to compare our results and the dif- 
ferences between the SFR and BHAR in terms of the an- 
alytic model given by Hernquist & Springel (2003) for the 
cosmic star formation history. Using detailed analytical 
arguments based on the dependence of the SFR on cos- 
mological growth of structure and the physics of star for- 
mation, they derived a fitting function for the SFR history 
given by 



p*(z) = p*(0) 



X 



l + a ( X -l) 3 exp(/?x 7 / 4 )' 



(12) 



where \ = [H{z)/Hq\ 2 / 3 , and where a, (3 and p(0) are 
parameters. In fact, Hernquist & Springel (2003) showed 
that this form (also drawn as the dotted line in Fig. 9) 
provides a more accurate fit to the simulation results for 
the star formation rate than equation (11). 

In particular, Hernquist & Springel (2003) demonstrated 
that the low redshift dependence of the star formation rate, 
p*{z) oc _ff(z) 4 / 3 , is primarily caused by the declining effi- 
ciency of gas cooling, which itself is related to the expan- 
sion rate of the universe, as measured by the evolution of 
the Hubble constant, H{z). At high redshifts on the other 
hand, the evolution is driven by the gravitational growth of 
structure. It is therefore not too surprising that the BHAR 
closely matches the SFR at high redshifts, but shows a 
different behavior at low redshifts, since we assumed that 
black hole accretion is not linked in the same way as star 
formation to gas cooling, but instead to stellar feedback 
processes, most notably to galactic outflows. In the bot- 
tom panel of Figure 9 we also show a fit to the BHAR 
using a function similar to equation (12) where, in order 
to fit the steeper decline at intermediate to low redshifts 
we have assigned the numerator to be oc \ instead of cx x 2 
(and a = 0.013, (3 = 0.12 and pbh(0) = 2.2 x 10" 4 p*(0)). 
At intermediate to low redshift the BHAR therefore scales 
as pBH oc p*(z) x H(z) 4 / 3 oc H(z) 8 / 3 ; i.e. roughly as the 
square of the cooling rate. This new scaling may arise as a 
result of the fact that gas is depleted by the joint effects of 
star formation and feedback. The BHAR then may scale 
in proportion to both the star formation (which scales as 
H(z) 4 / 3 in this regime) and the feedback processes (which 
themselves scale as the SFR), hence as the square of the 
cooling rate. In future work, we will investigate in more 
detail possible physical processes underlying this scaling 
by means of analytical arguments. 

6. CONCLUSIONS 

We have explored the hypothesis that black hole growth 
and activity at the centers of galaxies is primarily reg- 
ulated by the evolution of the total baryonic gas mass 
content of galaxies. We have used ACDM cosmological 
hydrodynamic simulations, which include a converged pre- 
scription for star formation, to follow the evolution of the 
baryonic mass in galactic potential wells. Using the nu- 
merical models, we have searched for correlations in the 
large-scale properties of galaxies that match those of the 
observed black holes. Note that with such an approach, 
resolving the properties of the actual small-scale accretion 
flows around black holes is not important. 

We have shown that the steep power law form of the ob- 
served Mbh — er correlation can be understood over an in- 
teresting regime of circular velocities and redshifts if there 



is a simple linear relation between the total gas mass (sub- 
ject to star formation and associated feedback) in galaxies 
and their black hole mass. We find that the observed black 
hole mass density is consistent with Mbh ~ 0.004/i _1 M gas ; 
i.e. the central black holes contain < 1% of baryons in 
galaxies. The total amount of gas in galaxies, and hence 
their black hole mass, may saturate simply in response 
to star formation, and in particular to supernova feed- 
back and galactic winds. By considering these as the 
only physical processes regulating the growth of black 
holes, we predict that the steep dependence of black hole 
mass on velocity dispersion is tight but not set in primor- 
dial structures, but only fully established at low redshifts, 
z < 2. Once the relationship is established, we find that 
the black hole mass is related to the dark matter mass by 
M BH /10 8 M Q ~O.7(M DM /1O 12 M ) 4 / 3 . 

We assume that all galaxies undergo a quasar phase with 
a typical lifetime, £q ~ a few xlO 7 yr (of the order of the 
Salpeter time) and show that the evolution of the quasar 
gas supply (taken as the ~ 1% fraction of the total gas) 
in spheroids is sufficient to explain, for the most part, 
the decrease of the bright quasar population at redshift 
z < 3. However, at redshift z <~ 1 the predicted shape of 
the luminosity function differs from observations implying 
that mass dependent variations in the typical accretion 
timescale or radiative efficiency are likely to be impor- 
tant. An additional possibility is that there is significant 
intrinsic obscuration (as suggested by recent Chandra con- 
straints on the hard XRB; e.g., Rosati et al. 2002) which 
may also be mass dependent. 

The quasar number density at high redshift in our model 
does decline, but under the simplest assumption of a red- 
shift independent t q , it can still exceed the SDSS space 
density of bright quasars at 4 < z < 6 (e.g., Fan et al. 
2001a, b) by a factor ~ 2 or more. However the observed 
i?-band number densities at z > 4 may be significantly 
biased low by dust along the line of sight (Pei 1995; Fan et 
al. 2001a), as is suggested by observations of high redshift 
quasars which are found to reside in dusty environments 
(e.g., Omont et al. 1996; Carilli et al. 2000, 2001). These 
biases in in B-band selection can be circumvented by a 
comparison of our predictions with the newly constrained 
comoving number density of X-ray selected quasars (Cowie 
et al. 2003). We have done this and find that the model 
yields the observed evolution of number density of X-ray 
faint quasars over the whole redshift range 1 < z < 6 al- 
though the statistical significance of these observations is 
still not very high. 

We further note that any freedom we have to change the 
luminosity function of quasars in our model comes down 
to our ability to select the accretion timescale £q, which 
is our only free parameter in the model. This is a fun- 
damental parameter which governs black hole and quasar 
evolution, and observational constraints on its value will 
make our model predictions more firm and comparisons 
with data even more interesting. In the next few years 
constraints on quasar lifetimes are likely to become avail- 
able from study of the clustering and the proximity effects 
of quasars mostly from the SDSS. 

It has been shown that models based on the Press- 
Schechter formalism, where the quasar luminosity func- 
tion is related to the dark matter halo mass, and quasar 
emission is assumed to be triggered by galaxy mergers, 
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can reproduce the statistical properties of the observed 
SDSS quasar luminosity functions if appropriate assump- 
tions for the quasar duty cycles are made (e.g., Wyithe & 
Loeb 2002; Volonteri et al. 2002; Hatziminaoglou et al. 
2002). However, in these models, the decrease in merging 
rates is not sufficient to reproduce the luminosity func- 
tions at low rcdshifts and additional modeling of the evo- 
lution of the gas fraction is then required (Kauffmann & 
Haehnelt 2000). This highlights the importance of prop- 
erly including the physical processes that govern the evo- 
lution of baryons trapped in halos as well as the growth of 
galaxies and their interactions. In this work we did not as- 
sume that quasar emission is triggered by galaxy mergers, 
but simply assumed that all galaxies that are star-forming* 
have an equal probability of hosting a quasar. Because the 
quasar lifetime tQ is short, there can be many generations 
of quasars in our model, but the number density of quasar 
hosts is significantly larger at high redshifts, where tQ be- 
comes a more significant fraction of the age of the universe 
than at low redshifts. 

We have used our model for the relation between gas 
mass and black hole mass to predict the evolution of 
the black hole mass density in galaxy centers. We have 
shown that if black hole mass is assembled by gas accre- 
tion, our predicted value for the black hole mass density 
Pbh ~ 2 - 4 x 10 5 M Q Mpc~ 3 (extrapolated to z = 0) is 
consistent with observational estimates for the local value 
(Yu & Tremaine 2002). We find that the large majority 
of black hole mass is assembled up to and preceding the 
peak of the bright quasar phase at z ~ 3 and that almost 
no further growth takes place at lower redshifts. 

We have also derived the evolution of the black hole ac- 
cretion rate density (BHAR) and showed that the major- 
ity of black hole accretion should occur in the high density 
environments at z <~ 4 — 5, in rough correspondence with 
the peak of the star formation rate (SFR) history in the 
simulations. This is consistent with the derived black hole 
mass density evolution and implies that the peak of the 
optically bright quasar phase occurs only when the largest 
black holes are already assembled. At very high redshift, 
the BHAR and the SFR evolve similarly, both primarily 
driven by the rapid gravitational growth of structure, while 
at low redshift, the BHAR declines much more rapidly. Al- 
though the peak of the BHAR occurs close to that of the 
star formation rate, its normalization is a factor of a few 
hundred lower than that of the cosmic SFR. It hence seems 
unlikely that black hole accretion plays a crucial role for 
the gas dynamics in galaxies, and it may even be relatively 
unimportant for the general process of galaxy formation, 
unless there is strong energetic feedback by active QSOs 
that affects galaxies in a significant way. In future work we 
will investigate the predictions of this model for very high 
rcdshifts and in particular its implication for the nature of 
the ionizing background. 

This work was supported in part by NSF grants ACI 96- 
19019, AST 98-02568, AST 99-00877, and AST 00-71019. 
The simulations were performed at the Center for Paral- 
lel Astrophysical Computing at the Harvard-Smithsonian 
Center for Astrophysics. 

*Note this choice does not affect our results. It simply excludes 
small-mass objects where star formation is likely not to be resolved. 
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